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Bosons and fermions, in the presence of frustration or background gauge fields, can form many- 
body ground states that support equihbrium charge or spin currents. Motivated by the experimental 
creation of frustration or synthetic gauge fields in ultracold atomic systems, we propose a general 
scheme by which making a sudden anisotropic quench of the atom tunneling across the lattice and 
tracking the ensuing density modulations provides a powerful and gauge invariant route to probing 
diverse equilibrium current patterns. Using illustrative examples of trapped superfluid Bose and 
normal Fermi systems in the presence of artificial magnetic fluxes on square lattices, and frustrated 
bosons in a triangular lattice, we show that this scheme to probe equilibrium hulk current order 
works independent of particle statistics. We also show that such quenches can detect chiral edge 
modes in gapped topological states, such as quantum Hall or quantum spin Hall insulators. 



The physics of fermions or bosons moving in back- 
ground gauge fields is of great interest in various con- 
densed matter systems such as quantum Hall liquids 
topological insulators [2 , quantum spin liquids [3 , and 
the cuprate superconductors [4]. Such abelian or non- 
abelian gauge fields, imposed externally or generated by 
strong correlation effects, can result in equilibrium charge 
or spin currents of electrons. For instance, a type-H su- 
perconductor in a magnetic field forms an Abrikosov vor- 
tex lattice that supports a periodic bulk current pattern 
formed by Cooper pairs swirling around each vortex [5^. 
A uniform magnetic field for lattice electrons can lead 
to topologically nontrivial states with a quantized Hall 
conductance and chiral edge currents 0. Such electronic 
charge currents in a solid produce their own character- 
istic magnetic fields, and can thus be probed by using 
magnetic microscopy or neutron scattering. These tools 
have been used to study vortices in type-H superconduc- 
tors [7 , to search for complex current patterns in the high 
temperature cuprate superconductors [8 , or to look for 
edge currents in purported chiral superconductors such 
as SrRu204 pi. Electronic spin currents in solids, by 
contrast, are harder to measure. A direct observation of 
the spin Hall effect in semiconductors involves driving a 
charge transport current and optically detecting the spin 
accumulation at the transverse edges of the sample [10]. 
Observing equilibrium spin currents is a more difficult 
challenge; only recently have experiments shown that the 
quantum spin Hall edge modes in two-dimensional HgTe 
quantum wells carry spin polarization [11 . 

Over the past few years, experiments in the field of ul- 
tracold atomic gases have also begun to study the effects 
of "artificial" orbital magnetic fields p!2HT4] and spin- 
orbit coupling [15 in the hope of creating new states of 
atomic matter. These experiments can potentially realize 
various topological phase transitions and a wide variety of 
states with equilibrium mass currents [16l[T7]. Such mass 



currents also arise in the presence of 'lattice shaking' 
[iHl |19] , the combination of Raman lasers and radio fre- 
quency fields [20^, or from populating higher optical lat- 
tice bands with bosons [211 122] , both of which lead to ki- 
netic frustration and possible spontaneous time-reversal 
symmetry broken superfiuids [23l [24]. Two-component 
bosons, in the presence of spin-orbit coupling and strong 
correlations, have recently been proposed to support 
Mott insulator states with complex magnetic textures, 
such as vortex crystals and skyrmion lattices [25H27] . 
Upon decreasing the Hubbard repulsion, such Mott insu- 
lators transition into superfiuids, which retain the mag- 
netic textures, with the magnetic order imprinting non- 
trivial Berry phases on the bosons and leading to intricate 
superfluid current patterns [25 . Spinless fermions with 
longer range repulsive interactions and frustrated hop- 
ping on the triangular lattice have also been recently pro- 
posed to realize states with spontaneously broken time- 
reversal symmetry and loop currents [28] . 

But, how can we experimentally deduce the equilibrium 
mass current patterns for such neutral atomic gases? 
This is rapidly becoming an important issue since cold 
atomic gases are poised to create a number of interesting 
condensed matter states using such gauge fields. Ex- 
periments on bosonic atoms use peaks in the boson mo- 
mentum distribution to infer the location of the boson 
dispersion minima induced by the presence of synthetic 
magnetic fluxes [T3^; it would thus be extremely valu- 
able to have a complementary technique that directly 
probes gauge invariant equilibrium mass currents induced 
by such background synthetic gauge fields or frustration. 

In this paper, we argue that the study of density dy- 
namics triggered by specific quantum quenches provides 
a powerful route to probing atom mass currents, and we 
present an extended discussion of this idea going well be- 
yond our previous work [29]. Our proposal to measure 
equilibrium atom currents induced by the presence of a 
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gauge field relies on measurements of the atom density, 
and is inspired by the significant experimental progress in 
measuring even lattice scale density modulations. Many 
such density mapping tools have been experimentally 
demonstrated in recent years, such as noise correlations 
[3T] , Bragg scattering [32 , which is analogous to X- 
ray scattering used to deduce the crystal structure of 
solids, and in situ microscopy which is similar to 
scanning tunneling microscopy at a crystal surface in the 
sense that both probe real-space lattice scale physics. 

One key idea we use is to make a specific quantum 
quench of the Hamiltonian that violates the steady state 
divergence- free condition on equilibrium currents [29|. 
This causes an imbalance between currents entering and 
leaving different sites of the lattice. As dictated by the 
continuity equation, this leads to characteristic density 
build-up or depletion with a specific pattern across the 
lattice, which reflects the initial currents in equilibrium. 
A "quasi-local" current probe of this type has been used 
in a recent experimental study of nonequilibrium dynam- 
ics in a (one dimensional) ID Bose gas [3T. In other 
cases, one can design suitable quenches that lead to spon- 
taneous macroscopic dipolar density oscillations, corre- 
sponding to center of mass oscillations of the atom cloud 
in the harmonic trap. In either case, imaging the subse- 
quent density variation across the lattice yields real space 
information about the initial currents. Thus, just as the 
usual time of flight images probe momentum informa- 
tion by studying real space atom positions after a time 
delay following release from a trap, our proposed scheme 
yields atom current information by converting them into 
density images after a time delay. Such quenches along 
with the underlying current patterns are schematically 
depicted in Fig. [l] for some of the examples explored in 
this paper. 

Our proposed scheme has the additional advantage of 
being independent of particle statistics, and is applica- 
ble to both bosons and fermions, as we illustrate here for 
both Bose superfluids and degenerate Fermi gases. More- 
over, it can also be used to probe both bulk currents and 
edge currents in the system. In addition to the systems 
explored here, we also expect that such quenches could 
also probe the current pattern in the recently studied 
chiral Bose Mott insulator ^35ll36], and, more generally, 
the dynamics can also be used to study spin currents of 
atomic matter, since one can experimentally probe the 
spin-resolved density in the lattice as demonstrated in 
recent experiments [37] , 

Finally, quantum quenches have long been of great in- 
terest in the context of such diverse and important issues 
as the approach to equilibrium in closed quantum sys- 
tems, defect production induced by tuning the Hamilto- 
nian across various quantum phase transitions, and ex- 
tensions of scaling and renormalization group ideas to 
dynamics across quantum phase transitions such as in 
the context of the Kibble-Zurek problem [38H45] . Our 
work, thus, additionally serves to bring together these 
two threads of research — synthetic gauge fields and 




FIG. 1: Illustration of the investigated scenarios: (a) Uni- 
directional quenches of the tunnel coupling in a square lat- 
tice with staggered (checkerboard) magnetic flux, (b) uni- 
directional quenches in a square lattice with a striped flux 
pattern and (c) bidirectional quenches in a triangular lattice 
with frustrated hopping. Strong and weak tunnel couplings 
are illustrated by thick and thin lines and the magnetic unit 
cell is highlighted. The quench dynamics for all three cases 
was investigated considering a Bose superfluid; case (a) was 
additionally studied for noninteracting fermions. 



quantum quenches — by suggesting that studying quan- 
tum quenches and quench-induced dynamics in the pres- 
ence of gauge fields would be a useful direction to pursue 
given the recent experimental and theoretical advances 
in these areas. Indeed, there have even been theoretical 
proposals to produce dynamical gauge fields in ultracold 
atomic systems [46^ 47^ , and a recent suggestion that one 
could use quenches in Bose-Fermi mixtures to simulate 
'string breaking' dynamics in a model of fermions coupled 
to such dynamical gauge fields [48 . 

A brief version of some of the results in this paper is 
contained in Ref. [29j. Here, we outline a general scheme 
to extract currents, expand on the various analytical re- 
sults for noninteracting cases and for interacting Bose 
systems, give further details on quenches for quantum 
Hall states, study a triangular lattice frustrated super- 
fluid, and explore the effect of random phase fluctuations 
imprinted on the initial state prior to the quench in or- 
der to show that weak thermal fluctuations in a Bose 
superfluid do not affect our central conclusions. 



I. GENERAL SCHEME AND MODELS 

Consider a general (i-dimensional lattice system that 
carries local currents j(r) in equilibrium. In order to 
uncover a specific component of the current, say jx = 
j • X, let us make a quench of the Hamiltonian at time 
t = that instantaneously turns off the hopping along 
all transverse directions. The currents j(r) remain un- 
changed at the instant of the quench; the quench itself 
induces no extra currents. However, at all subsequent 
times, the density and current will evolve in an effec- 
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lively one-dimensional system along the x direction. The 
continuity equation following the quench takes the form 

-^^=J.(r + x/2)-j,(r-x/2), (1) 

where the right-hand side denotes the lattice divergence 
of the current. 



A. Short time analysis 

If we focus on short times after the quench, the change 
in density (5n(r, t) will be dictated by the initial equi- 
librium currents in the unquenched direction, indepen- 
dent of the details of the Hamiltonian or the statistics 
of the particles making up the fluid. Therefore, at short 
time e after the quench and in Fourier space, we can set 
An(q, e) ^ -2ze sin(g'a,/2)ja,(q, 0), with ja,(q,0) being the 
Fourier transform of the pre-quench equilibrium current. 
For this short time analysis to be valid, we must choose 
1/e to be comparable to the tunneling rate, but much 
shorter than the frequency of the subsequent density os- 
cillations discussed below. 

For Qx 7^ 0, we can invert this relation to obtain 

J.(q,0) - .] /^J n{q,e). (2) 

2iesm[qx/2) 

Thus, a measurement of the excess density Sn, which 
builds up shortly after the quench and decomposing it 
into its spatial Fourier components, corresponds to a de- 
termination of all the nonzero Fourier components of jx- 

To determine the Fourier components of the density, 
one could resort to tools such as Bragg scattering [32 , 
noise correlation measurements [30l ISTJ or superlat- 
tice aided band-mapping techniques [50] . Assuming sim- 
ple current patterns, only a few Fourier components will 
be nonzero. For low filling < 1, in situ imaging after 
freezing out all atom hopping provides the most direct 
measurement [33l |5T^. For larger filling, however, ad- 
ditional effort has to be made to overcome the parity 
mapping inherent to this method. 

For Qx = 0, the above inversion fails. This component, 
however, corresponds to a uniform current offset and can 
be detected in the presence of an external trapping po- 
tential by monitoring the change in the center-of-mass 
position after short times e. 

Aside from the short-time dynamics following a com- 
plete quench, it is also useful to study the long-time dy- 
namics and the dynamics following a weak quench; both 
these issues are amenable to analysis and experiments, 
and are interesting in their own right. At the very least, 
such an analysis of the dynamics allows us to address the 
issue of how long the system needs to evolve before the 
measurement to ensure there is an experimentally mea- 
surable density accumulation. In order to make progress 
on this front, we focus on specific model Hamiltonians. 



B. Models 

We are interested in applying the general scheme out- 
lined above to uncover the underlying current patterns 
of interesting and experimentally relevant examples of 
many body states of bosons and fermions. We therefore 
focus on the following models of Bose superfluids: {i) 
bosons on a square lattice in a staggered, checkerboard- 
like magnetic flux pattern as shown in Fig. [ija), (n) 
bosons on a square lattice with a striped magnetic flux 
pattern as realized in Ref. 13 (Fig. [ijb)), and (m) 
a triangular lattice model of frustrated bosons mov- 
ing in a staggered flux pattern as realized in Ref. [18] 
(Fig. [ijc)). In addition, we study two models of non- 
interacting fermions: {iv) non-interacting fermions on a 
square lattice in a staggered magnetic flux (Fig. [ija)), 
and (v) the Hofstadter model of fermions with a uni- 
form flux on a square lattice, which results in an integer 
quantum Hall phase. As depicted in Fig.[l} we study sud- 
den quenches where we turn off (or weaken) the hopping 
along all directions except one. For the cases (i), (ii) 
and {iv)^ the quench we study corresponds to turning off 
or weakening the hopping along one direction, leading to 
sublattice-density oscillations. In case (ni), the quench 
turns off the hopping along two of the three bond direc- 
tions, resulting in macroscopic dipole oscillations of the 
atom cloud in a trap. For the quantum Hall case (v), 
a unidirectional quench is shown to lead to quadrupole 
oscillations dominated by chiral edge currents. 

In order to model the density dynamics of these sys- 
tems, we resort to two approaches. For bosons, we study 
examples with repulsive contact interactions modelled by 
a Hubbard Hamiltonian of the schematic form 

HbH = -Y1 JrAK'+Jl^r b%+J E ^rblbA- (3) 
r,r' r r 

Here, the complex hopping amplitudes (transfer inte- 
grals) Jj. encode the artificial fluxes, U is the on-site 
Hubbard repulsion, and = Vb(x^ + ?/^) is a harmonic 
trap potential. We analyze these boson models using a 
Gross-Pitaevskii (GP) approach, which replaces by a 
"condensate wavefunction" ^j.. This leads to an equilib- 
rium energy functional 

r,r' r r 

which must be minimized to obtain the initial equilibrium 
superfluid ground state. Starting from this ground state 
of the pre-quench Hamiltonian, we suddenly decrease the 
tunneling amplitudes from their initial values, (J^, J^), 
to their final values, (J/, J^), at time t = and study 
the subsequent time evolution of this state. The quench- 
induced dynamics is then obtained by solving the time 
dependent GP equation 
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Henceforth, we set h = 1. While we present an analyti- 
cal discussion of the quench- induced dynamics in uniform 
systems, we also present numerical solutions for Bose 
superfluids in a harmonic trap. Specifically, the equi- 
librium state is obtained by numerically minimizing the 
GP energy functional, while the post-quench dynamics 
is obtained by numerically solving the time-dependent 
GP equation. Details of the numerical procedures are 
contained in Appendices A and B. 

For fermions, we restrict ourselves to noninteracting 
(spinless) examples for which we can diagonalize the 
Hamiltonian either analytically or numerically for large 
systems. These fermion Hamiltonians schematically take 
the form 

^ = - E •^r.r' /rVr' + E /rVr ■ (6) 

r,r' r 

We again imagine quenching the fermion hopping at time 
t = 0, with the time evolution being governed by the 
time-dependent Schrodinger equation. Knowing all the 
pre-quench and post-quench eigenstates and eigenvalues 
is then sufficient to reconstruct the dynamics of various 
observables after the quench. 



II. BOSE SUPERFLUID IN A STAGGERED 
MAGNETIC FLUX BACKGROUND 



Let us restrict ourselves to fiux values < <j ) <tt. W e find 
that the minimum eigenvalue, Ak = ~\/^k + 7k' ^^^^ 
occurs at k = (0,0). For = Jy = J in equilibrium, 
this is given by 

Ao = -4Jcos^. (10) 

The GP wavefunction in the absence of a trap is given by 
the wavefunction corresponding to this minimum eigen- 
value, = \^{uo + ivoTjr)^ whcrc 

1 / r \^/^ 

- = ('-£)"■ <'^> 

and r]r = (—1)^^^. In this initial equilibrium state, the 
density |^rP = is uniform and there is an alternating 
checkerboard pattern of circulating currents on the ele- 
mentary square plaquettes. The magnitude of this stag- 
gered current is given by 4JiioVo^o on each bond. 

Since the density in this state is uniform, this wave- 
function continues to be the ground state of the full GP 
equation in the absence of a trap. In later subsections 
where we present a numerical solution to the GP equa- 
tion in the presence of a trap, the currents and densities 
are nonuniform. 



We begin by studying a weakly interacting superfiuid 
of bosons on a 2D square lattice, described by the Bose- 
Hubbard model Eq. (|3|. We take Jry = J*' ,r c>nly 
for nearest neighbors, and choose Jr,r+x = Jx i"eal and 
Jr,r-\-y = Jy exp[z(— 1)^+^^/2] . This yiclds staggered mag- 
netic fiuxes, that pierce the elementary square pla- 
quettes in a checkerboard pattern [see Fig.jlja)]; a route 
to realizing such a fiux pattern has been proposed previ- 
ously [54^ . 



A. Equilibrium state in the absence of a trap 

For weak interaction, U ^Jx-,Jyi we solve for the equi- 
librium ground state by minimizing the GP energy func- 
tional Eq. Q [52 . In the absence of a trap {Vr = 0) 
we first diagonalize the kinetic energy in the Hamilto- 
nian ([3|. In momentum space, the kinetic energy takes 
the form 



H 



E 

l^eRBZ 



(7) 



where r^'^ are Pauli matrices, Q = (7r,7r), and RBZ 
denotes the reduced Brillouin zone due to the unit cell 
doubling resulting from the fiux. Here, we have defined 

= — 2 Jx cos kx — 2 Jy cos{(j)/ 2) cos ky (8) 
7k = — 27^^ sin((/)/2) cos /c^ . (9) 



B. Exact analysis of a quench for noninteracting 
bosons with no trapping potential 

For noninteracting bosons, it is simple to analyze the 
quench dynamics in the absence of a trap, since we ex- 
plicitly know the energies and eigenstates before and af- 
ter the quench. Specifically, let the equilibrium time- 
independent wavefunctions in the pre-quench and post- 
quench Hamiltonians be given by ^/no{uo + ivoVr) and 
^/no{uo + ivoVr) respectively, where the coefficients of 
the uniform and staggered components are determined 
from Eq. (11) and Eq. (12). We can then write the post- 



quench time-dependent wavefunction in the form 



^(r,t) 



+ (3{vo - iuoVr)e'~^''] , (13) 



where Aq is the lowest energy eigenvalue of the post- 
quench Hamiltonian, and 



a 



{uqUo + VqVo) 
{uoVo - VqUo). 



This leads to a time dependent density 

n{r,t) = \^{r,t)f 

= no 1 — 2a/37^r sin(2|Ao|t) 



(14) 
(15) 



(16) 
(17) 
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which exhibits staggered modulations on top of the uni- 
form background average with a frequency 2|Ao| and an 
amphtude, which depends on the degree of the quench. 
For a weak quench, where uq^vq are close to uq^vq^ the 
amplitude is small; however the amplitude can be signif- 
icant for a strong quench. 



C. Approximate analysis for interacting bosons in 
the absence of a trap 

Let us now consider the effect of interactions on a 
weak quench in the absence of a trap, where sud- 
denly decreases from J to J -\- 5J < J at time t = 
while keeping Jy = J constant. (Note that since we 
start from the isotropic case, and since we are studying 
currents and densities, which are both gauge invariant 
quantities, we would get exactly the same results for a 
quench along the ^/-direction.) Because the quench con- 
serves crystal momentum in the reduced Brillouin zone, 
we can write the post-quench wavefunction in the form 
^(r,t) = A{t)^B{t)r]r, where A{t) and B{t) denote time 
dependent complex coefficients of the uniform and stag- 
gered component of the wavefunction in real space. The 
full time-dependent GPE then reduces to a pair of non- 
linear ordinary differential equations for A{t) and B{t)^ 
given by 

HA 

'Itt ^ " '^'^ ^ U{A\A\^ + 2A\B\^ + 52 A*), (18) 

= ^70^ - ^0^ + U{B\B\^ + 2B\A\^ + A25*),(19) 

where 70 = — 2Jsin(0/2) and eo = — 2(J + 5 J) — 
2Jcos(^/2). It is easy to check that the total density 
does not change, since d/dt{\A\^ + I^P) = 0- However, 
as discussed in detail in Appendix C, the staggered dif- 
ferential density AnAsit) = 2 (A* 5 + 5* A) can be shown 
to approximately obey the simple harmonic equation 



d^AuAB 
dt^ 



(20) 



where = 4[Aq + /7no(eoeo +7o7o)/|Ao|]. Using the ini- 
tial condition on the equilibrium currents, the solution 
to this can be written in the intuitive form 

AnAB{t) = 4 (^^^^sinm (21) 

where X is the initial current on each bond, given by 
4J^no sin((/)/2) 



1 = 



|Ao| 



(22) 



SJ/J represents the fractional change in the hopping 
along the x-direction (which is also the instantaneous 
fractional change in the current along the x-direction in- 
duced by the quench), and the factor-of-four in the front 
arises from two x-bonds having been weakened by the 



quench. These sublattice density oscillations thus di- 
rectly reflect the presence of staggered currents in the 
initial state; the amplitude of these oscillations depends 
linearly on SJ for a weak quench, while the frequency 
increases with increasing interaction strength. Knowing 
J and (5 J, a measurement of AnAB{t) and its oscillation 
frequency Q would thus provide quantitative information 
about the initial equilibrium current X, which can be 
compared with the theoretically expected value quoted 
above. 

In case experimental imperfections are too strong to 
observe oscillations with a well-resolved frequency, the 
initial build-up of the density pattern might be used to 
extract X. While this approach is unaffected by e.g. spa- 
tial variations of 1^, it rests on the ability to detect small 
changes in the sublattice populations. 



D. Numerical study of quench dynamics in the 
presence of a trap 

Having understood the underlying quench dynamics of 
the staggered flux state in the bulk, we now reintroduce a 
harmonic trap potential. The equilibrium state is solved 
self-consistently and leads to a superfluid ground state 
with staggered loop currents as shown in Fig. [2|a) for a 
system with linear length L = 22, Vq = 0.07J, U = 0.2 J 
and an average filling factor of no = 4. The smooth den- 
sity profile of the ground state reflects the trap potential, 
but it does not reveal the currents induced by the gauge 
field. 

Starting with the equilibrium ground state, we perform 
a quench along the the x-direction and study the subse- 
quent density dynamics. As directly seen in Fig. |2|c,d), 
the condensate develops striking checkerboard oscilla- 
tions at t > that reflect the underlying current or- 
der. These oscillations can be monitored by the con- 
trast of the spatial sublattice density modulations Cab = 
[NA{t) - NB{t)]/N shown in Fig. ^h) (with TV being 
the total number of bosons) . Information about the di- 
rection of circulation on a plaquette is easily discerned 
from the density pattern established after a short time 
period; since the quench is in J^^ the initial build up 
of density is on sites that have currents flowing into 
them along the strong J-bonds oriented along the y- 
direction. After a short time has passed, the density 
buildup reaches a maximum and the flow is reversed, re- 
sulting in "plasma oscillations" between the two checker- 
board patterns. The frequency of these oscillations scales 



as ~ (Aq + Uno{eoeo + 7o7o)/|Ao|)) as shown earlier; it 
thus varies slowly with position due to the inhomogeneity 
of the density in the trap. 

In order to compare the numerically computed dy- 
namics with the analytical results of the previous sec- 
tion, we compute the local sublattice density contrast 
C^^^ = {nA(t) — TiB(t))c', over a region of 4 x 4 sites in the 
center of the trap. We find that C^ab^^) oscillates with a 
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FIG. 2: (Color online) Density pattern of 2D Bose super- 
fluid in a staggered flux gauge field following a quench, (a). 
Initial density profile and (inset) current pattern of conden- 
sate ground state at half- filling, (b). Time dependence of 
the sublattice density contrast CAs(t) for a staggered flux of 
101 = 7v/2 and U = 0.2 J for different cases: (i) Trap average 
of Cab (t) following a complete quench from J* = J to = 
(top, red), (ii) Average of CAB{t) over a central 4x4 region 
(centre, green) following a quench from J* = J to = 0, and 
(iii) Trap average of CAsit) following a partial quench from 
Jx — J to JI — 0.25 (bottom, blue), (c - f). Change in lo- 
cal density, 5n (relative to original density), at different times 
following a quench = J ^ = OJ with U = 0.2J. The 
marked region in (c) and (d) indicates the central four-site 
plaquette. 



significant amplitude (up to ~ 25% contrast) and it can 
be fitted with the form in Eq. ( [2T| . The value of the cur- 
rent extracted from such a fit is 14.8 J, which is remark- 
ably close to the value obtained from the analytical ex- 
pression 4nJcos((/)/8) sin((/)/8) c:^ 14. 9J (the parameters 

are the central density u'q^ = 19.4 and = 7r/2). The ex- 
tracted value of the oscillation frequency, Q = 6.72J~^, 
also agrees well with the analytic result of 6.68 given 



by Eq. (20). It should be noted that when the system is 
taken altogether, the early-time dynamics would provide 
a better fit to the average current, which is dominated 
by the central region. 

The sublattice density difference when integrated over 
the entire trap exhibits some degree of dephasing and 
damping due to the density inhomogeneity; nevertheless. 



given our finite system size, the sublattice density oscil- 
lations persist out to fairly long times, t J ^ 1, as seen 
in Fig. [2|b). At these times, we find additional long- 
wavelength modulations superimposed on the checker- 
board density pattern (see Fig.|2|e,f)). Prominent spher- 
ical density waves emanate periodically outward from the 
centre, which we attribute to the spatial variation of the 
"plasma frequency" resulting from the radial variation 
of the density |^rP in the trap. Furthermore, the cloud 
shape shows oscillatory distortions into an ellipse due to 
the anisotropy of the final tunneling < Jy. 



E. Effect of random noise in the initial state 

As a simple test of the robustness of the quench pro- 
cedure to condensate depletion, we compute the quench 
dynamics of a state with imposed random phase fluctu- 
ations so as to mimic thermal fluctuation effects. We 
first add a small random and uncorrelated phase shift to 
each site, which is chosen in the range (0,^^max)- This 
'random' state is evolved for a long period of time ac- 
cording to the (unquenched) GP equation to allow the 
system to equilibrate into a viable 'thermal state', which 
now supports correlated phase and density fluctuations. 
Starting from this thermal state, we next perform the 
sudden quench by evolving this state according to the 
quenched GP equation, as before, and analyze its dy- 
namics. We consider two cases, one with small fluctua- 
tions (5^mlx = 0.5rad and one with moderate fluctuations 
50m.L^ = l.Orad. To the extent that these fluctuations 
lead to states that mimic typical states from a thermal 
ensemble, both realizations of phase fluctuations result 
in superfluid states well below the Berezinskii-Kosterlitz- 
Thouless transition; details are given in Appendix D. 

As shown in Fig. [sj^a), the initial sublattice density 
contrast in a given state with noise, or averaged over 
several such realizations, is similar to that of the clean 
system; however, whereas the oscillations persist for a 
long time in the ground state, the oscillations in such 
a'thermal state becomes incoherent after a few periods. 
Furthermore, even at short timescales, we find that: (i) 
the amplitudes are no longer equal to that of the ground 
state quench, and (ii) the oscillation frequency is slightly 
shifted compared to its zero temperature result. 

A further impact of thermal fluctuations is found in the 
broadening of the spectral peak in the structure factor 
shown in Fig. IsFb-d). (Note, that the strong momentum 
peaks about (0,0) have been removed for clarity.) The 
large single spectral peak at (tt, tt) in the clean system, 
shown in Fig. [sjb), is replaced by a broader peak in the 
noisy systems shown in Fig.[3](c,d), where SOmax is equal 
to 0.5 rad and 1.0 rad in (c) and (d), respectively. Despite 
the fluctuations and broadening, the (tt, tt) peak at early 
times can still easily be discerned above the background. 
We observed that the real-space density pattern also dis- 
plays a discernible checkerboard-like tendency even in the 
presence of moderate thermal noise. 
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FIG. 3: (Color online) Dynamical density pattern of interact- 
ing bosons on a 2D square lattice in a staggered flux gauge 
field following a quench in the presence of noise, (a) Compar- 
ison of the time dependence of the sublattice density contrast 
CAB{t) for a staggered flux of = 7v/2 following a quench 
of J* = J to = without noise (solid red) and with 
noise (dashed blue). Noise is incorporated by including a 
random local phase fluctuation in the range {0,60max) with 
SOmax = 0.5 at each site, and evolving the system according 
the CPE for a long period of time before implementing the 
quench. Normalized structure factor at tJ = 0.25 for a sys- 
tem (b) without noise and (c,d) with the presence of noise. 
Noise is incorporated by adding a random phase, with fluctu- 
ation magnitude (c) SOm\x — 0.5 and (d) SOm\x — 1.0, locally 
to each site and allowing the system to equilibrate. Note, 
the momentum peak close to (0, 0) associated with the aver- 
age density distribution has been removed for clarity and the 
axes in (b - d) are equal. 



III. BOSE SUPERFLUID IN A STRIPE 
SYNTHETIC FLUX BACKGROUND 

We next consider the Bose-Hubbard model in the pres- 
ence of a striped magnetic flux pattern as realized in 
Ref. [13 (see Fig. [lib)). We choose Jr,r+y = Jy and 
Jy,y+x = exp[z( — lj^(/)y], so that we enclose fluxes 
through each plaquette that lies along a stripe in the y- 
direction, and solve for the equilibrium ground state by 
minimizing the GP energy functional for Jx = Jy = J ^ 
and for weak interactions U = 0.2J, with L = 22, 
Vb = 0.07 J, and an average filling no = 4. We find a 
superfluid with vertically striped loop currents depicted 
in Fig. |4]^a), which resembles a stripe pattern of 'long 
vortices' that are highly elongated along the ^/-direction. 
Again, the smooth equilibrium density pattern is reflec- 
tive of the underlying trap potential but it reveals no 
information about the underlying currents. 

Upon quenching J^^ the superfluid generates a density 
pattern that strikingly reflects the underlying equilibrium 
striped currents. Each vertically elongated loop forms 
four quadrants of alternating high and low density, giv- 
ing rise to an oscillatory quadrupole moments as can be 
seen in Fig. |4]^b-d). Evidence for striped density pattern 
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FIG. 4: (Color online) Dynamical density pattern of interact- 
ing bosons on a 2D square lattice in a stripe-like magnetic flux 
pattern following a quench, (a) Initial density profile and (in- 
set) current pattern of condensate ground state at half-filling. 
( b - d): Change in local density, n, at different times follow- 
ing a quench J* = J ^ = with [/ = 0.2 J. The circled 
region indicates the central elongated vortex with a nonzero 
quadrupole moment, (e) Density structure factor normalized 
to the peak maximum and with the momentum peak at (0, 0) 
associated with the average density removed for clarity, (f) 
Aspect ratio of the cloud as function of time for various fluxes. 



can also be found in the structure factor shown in Figure 
|4je). The small momentum modes about (0,0) that are 
attributed to the average density are subtracted off for 
clarity, but we emphasize that it is not necessary to know 
the original density distribution before the quench for any 
of the analysis. The structure factor shows two dominant 
spectral peaks in proximity to (tt, 0) but slightly shifted 
by q = (0,±^). This small momentum shift is a reflec- 
tion of the additional long wavelength component that 
originates from the anti-nodal line in the density pattern 
running along the centers of the elongated quadruples at 
y = 0. Since this q scales as 1/L, the two peaks will 
merge toward (tt, 0) for larger system sizes. 

Quenching Jy rather than leads to similar early- 
time density patterns; however, the oscillatory dynamics 
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that occur are much faster. In addition, the density os- 
cillations in this case rapidly decohere when compared 
with the quench in J^- These observations hint that a 
very different mechanism is at play for the two quenches. 

To better understand the dynamics, it useful to con- 
sider the infinite system without a trap. The initial den- 
sity pattern can be shown to be uniform with alternat- 
ing stripe currents fiowing up and down along the bonds 
that lie along the ^/-direction, which has a magnitude of 
2J^sin(^/2) for (j) < 7r/4. In this limit, neither quench 
(in Jx or Jy) leads to any density modulation and no 
information pertaining to the current can be obtained 
from such quenches. This situation occurs for any sys- 
tem with constant current along a given bond direction. 
Despite this, information about the current can still be 
inscribed onto the density when the quench is performed 
in the presence of a trap. 

Consider again the stripe flux state in the presence of 
a trap. Close to the center of the trap the current again 
alternates direction and flows along the Jy bonds, how- 
ever, there are now strong edge currents in the vicinity 
of the cloud boundary (in the above simulation, the edge 
currents die off as the center of the trap is approached, 
but are still present due to finite size affects). When the 
Jx hopping is quenched, the bosons now continue to fiow 
in the direction of current and, instead of leaking into the 
edge current, will fiow up the trap potential. This gener- 
ates density oscillations along the vertical chains whose 
frequency is determined by the trap profile and the effec- 
tive mass of the bosons. Since, the initial direction of the 
current fiow alternates along each of the vertical chains, 
the phase difference between the density oscillations of 
neighbouring chains is tt, thereby producing the stripe 
pattern. In contrast, a quench in Jy directly probes the 
edge currents in the system and, although the density 
pattern is again striped, the oscillation frequency is not 
determined by the trap potential. 

The least demanding experiment is a measurement of 
the aspect ratio of the cloud, ^D^'^ {t)/Dy2(t) where 
Dx'2 {y2){t) = '^j,n{T^t)x'^ (?/^), as a function of time (see 
Fig. pFf)). Notice that only the oscillation amplitude 
is discernibly affected by the value of the magnitude of 
the fiux per plaquette, while the oscillation frequency is 
essentially governed by the trap frequency and thus is 
practically independent of fiux. The variation of the am- 
plitude refiects the differences in the initial currents for 
different fiux values. 



IV. TRIANGULAR LATTICE FRUSTRATED 
BOSE SUPERFLUID 

A similar situation arises for a triangular lattice in the 
presence of a staggered fiux state (see Fig. [ijc)). Con- 
sider the initial state to be the ground state of a sys- 
tem with = 7r/2 fiux per plaquette and further Ji = 
J2 = Js = Ji where Ji, J2, and J3 are the magnitudes 
of the tunnel couplings along the (1,0), (1/2,y^3)/2) 
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FIG. 5: (Color online) Density dynamics for spinless 2D 
bosons in a staggered flux gauge field on a triangular lattice 
after simultaneously quenching the J2 and J3 bonds, (a) Ini- 
tial density pattern and (inset) current pattern at half filing 
with flux 101 = 7r/2. (b) The dipole moment as a function of 
time for various parameters. Distances are measured in units 
of the lattice constant. (c,d) Density proflle at the indicated 
times showing the displacement of the cloud. The crosshairs 
indicate the center of the trap. 



and (— 1/2, Y^3)/2) bond directions, respectively, and 
U = 0.2 J. As shown in the inset of Fig.jsja), the ground 
state currents flow along the bond directions, as they do 
for the inflnite system without a trap. Hence, we again 
rely on the combined effect of the quench with the trap- 
ping potential to transcribe information about current 
flow onto the density proflle. 

Since there are now three unique bond directions, it is 
necessary to quench the hopping along two of the bond 
directions simultaneously, which we choose to be J2 and 
J3. Here, the dynamics is very similar to that of the 
Jx quench for the stripe flux considered previously. Each 
chain decouples and the initial current causes the density 
to oscillate in the trap potential. This time, however, the 
direction of current flow is the same along neighbouring 
chains, resulting in a uniform oscillation of the entire 
cloud along the unquenched bond direction. (Note that 
unlike the checkerboard case, this current pattern is at 
q = 0, so that the quench does not produce any density 
modulations with nonzero Fourier components.) 

This oscillation is best observed by monitoring the time 
dependence of the dipole moment Dx = '^j,n{T^t)x^ as 
shown in Fig. jsjb); a larger equilibrium current in the 
initial state will lead to a larger amplitude for such center 
of mass oscillations following a quench. Notice, however, 
that the oscillation frequency is nearly independent of the 
magnitude of the quench. This substantiates the idea 
that the oscillation frequency set by the trap stiffness 
and is independent of the initial current. Furthermore, 
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the oscillation frequency is independent of the interaction 
strength, as we confirm numerically in Fig. jsjb), since 
it involves only center-of-mass oscillations. One should 
keep in mind that the oscillation frequency can shift if the 
initial flux per plaquette is changed. This is expected 
because the bosons will initially condense into a state 
with different crystal momentum and lattice effects will 
change the effective mass of the condensate. 



V. SPINLESS FERMIONS IN A STAGGERED 
MAGNETIC FLUX 




Motivated by our study of quench-induced density 
dynamics for bosons, we next turn to noninter acting 
fermions in a staggered flux background [53 . We study 
the Hamiltonian H^f = - X]r,r' Jr.v'fUv'^ where Jr.r+x = 
J and Jr,r-\-y = exp[z(— 1)^+^0/2], leading to stag- 
gered checkerboard fluxes ±0 [5T. To make analytical 
progress, we ignore the harmonic trap in the discussion 
below. As we have seen previously for bosons, and as dis- 
cussed below, the trap does not qualitatively affect our 
conclusions, and we can also directly apply our results 
to the central region of the trapped gas. In momentum 
space, the Hamiltonian takes the form 



i^sf = ^'l^k<(cos^kr" + sin^kr^)^^. 



(23) 



where Q = (7r,7r), = (/k^/k+g), and r^'^ are Pauh 
matrices. The prime on the momentum sum implies that 
only momenta in the reduced B rillouin z one are included. 
Here, we have defined l^k = V^k + 7k' cos^k = ^k/^k, 
and sin6>k = 7k/^k, with 



Sk = — 2(Jcos /Ca^ + cos — cos /c^) (24) 



7k = — 2 Jjy sin — cos ky. 



(25) 



This leads to mode energies ±l^k in the initial state. 

Imagine fermions initially filled into negative energy 
states — r^k up to a Fermi energy Ep^ and then quenching 
Jx from J* J I at time t = 0. Such a translationally 
invariant quench ensures that different momentum pairs 
(k, k + Q) stay decoupled from each other. Nevertheless, 
this quench instantaneously changes ^ S\^^ and 7k 
7k, so that we modify (l^k^^k) ^ (^k,^k)- 

This means that while the initial quasiparticle occu- 
pation numbers are set by the initial dispersions and the 
chemical potential, the subsequent dynamics is then de- 
termined by the final Hamiltonian. Since the final Hamil- 
tonian is also translationally invariant, the various mo- 
mentum states stay decoupled after the quench, but un- 
dergo the analogue of "spin precession" in the two-level 
(k, k + Q) space. 

To compute the density modulation between the two 
sublattices at a subsequent time, AnAsit) = {tia — ^b), 



FIG. 6: (Color online) Time dependence of the sublattice den- 
sity difference AnAs(t), for noninteracting spinless fermions 
on a 2D square lattice at a filling of no = 0.4, flux = 7r/2, 
following a quench from (upper panel) J* = J ^ = J or 
(lower panel) J* — J ^ J I = 0.5 J. 



we write 



M 
2 

M 



where M is the number of lattice sites. Carrying out the 
algebra, details of which are given in Appendix E, we find 

AriABit) = J^Y^ sin(^k - 4) sin(2(^kt) (27) 

{kjocc 

where the momentum sum runs only over initially occu- 
pied states in the reduced Brillouin zone. 

A numerical evaluation of the sum allows us to plot 
the sublattice density oscillations, shown in Fig. [6] for 
J* = J, JI = O.OJ and = 0.5J, a fermion density 
of no = 0.4 per site, and various staggered flux values. 
These oscillations exhibit multiple frequencies due to the 
large number of occupied fermion modes. However, over 
the entire range of displayed fluxes, and a wide range 
of densities no ^ 0.3 - 0.5 near half-filling, we find that 
the dominant oscillation frequency arises from initially 
occupied states near k=(0,7r) due to a van Hove singu- 
larity in the density of states. Picking this single mode 
k = (0, tt) in the above momentum sum leads to an esti- 
mated, nearly density-independent, dominant oscillation 

frequency (2^2*) ^ 4^J2 + (j|)2 - 2 Jj| cos f. Both, 
the flux dependence of this oscillation frequency for a 
partial quench, and its flux independence for a complete 
quench with = 0, are in quantitative agreement with 
the numerical data in Fig. |6] The larger density of states 
near k = (0, tt) also enhances the signal amplitude for fill- 
ings close to no = 1/2. The weak density dependence of 
AnAB{t) over a range of fillings indicates that trap in- 
duced inhomogeneities will not significantly affect these 
oscillations. 
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VI. TOPOLOGICAL STATES WITH EDGE 
CURRENTS 

Finally, we turn to gapped yet topologically nontriv- 
ial states such as quantum Hall insulators, Chern band 
insulators, or quantum spin Hall insulators, all of which 
have bulk gaps but support topologically protected edge 
currents. At a fundamental level, Chern band insulators 
are no different from integer quantum Hall states in the 
Hofstadter model, as discussed in recent work [55 . Both 
systems involve fluxes threading through plaquettes of 
the lattice with no net flux over an appropriately defined 
unit cell; in the Hofstadter model, this unit cell is the 
magnetic unit cell. 

Proposals to obtain such fluxes in experiments on cold 
atoms exist in the literature [34l [56H58] . The simplest 
models of 2D quantum spin Hall states or topological 
insulators, such as the Kane-Mele model [59 , may be 
viewed as two independent copies of quantum Hall insu- 
lators or Chern band insulators, with the two copies being 
labelled by a well-defined spin quantum number (equiv- 
alently 'hyperfine state' for an atom) and experiencing 
opposite magnetic fluxes. Much of the physics we discuss 
below, which involves studying the density dynamics fol- 
lowing a quantum quench, will then be applicable to such 
quantum spin Hall states if one can experimentally probe 
the density of each spin species. 

While recent work has focused on extracting the non- 
trivial band topology from time-of-flight measurements 
[60l[6l] or spectroscopy of the edge modes [62 , 63 , here 
we explore density dynamics induced by the unidirec- 
tional quench for lattice fermions in a uniform magnetic 
field. For concreteness and reasons of simplicity, we con- 
sider fermions on a 2D square lattice with a uniform 
magnetic flux (/) = 27r/3 per plaquette. Similar uniform 
flux configurations have recently been established in cold 
atom experiments by rotating the optical lattice [BH [65] . 
The resulting particle-hole symmetric Hofstadter spec- 
trum p6 has three non-overlapping bands, with Chern 
numbers +1, —2, +1, so that 'band insulators' with some 
bands being completely filled support a nonzero quan- 
tized Hall conductance, and chiral edge currents, yielding 
lattice versions of integer quantum Hall (QH) states in 
the continuum [6 . 

We begin by numerically diagonalizing the Hamilto- 
nian Hq^ = - Jrv'flfv' with Jr,r+£ = Je*^^ and 
J. 



r,r+y 



Jy^ for 



27r/3, with open boundary condi- 
tions on a L X L system, and fill up the lowest band (and 
some edge modes) to get a fermion filling no = 1/3. We 
find that the ground state bulk density is uniform (see 
Fig. [t] for t = 0) and supports edge currents confined to 
an "edge layer" of thickness ~ 2 — 3 lattice sites where 
the density also slightly deviates from its bulk value. We 
next track the density dynamics following a quench from 
Jl = J to JI < J, which is easy to study once we com- 
pute the initial and final spectrum and eigenstates. We 
note that the gauge choice for the magnetic field (i.e., how 
exactly to include the vector potential) is unimportant — 




FIG. 7: (Color online) Density dynamics for spinless fermions 
in the lowest Hofstadter band with flux = 27r/3 per plaque- 
tte on a square lattice following a quench from = J to 
Jx = 0.0 J. (a) The scaled quadrupole moment Qxy{t)/L^ = 
(1/L^) xyn(r, t) for various system sizes versus the scaled 
time tJ/L. (b - f) Stripe- like density modulations (for L = 24, 
plotted as n(r, t) — no) moving from the ^-edges into the ini- 
tially incompressible bulk at different indicated times. 



we could equally well have them along the x-bonds. 

Viewing the chiral edge currents as analogous to that 
arising from a 'vortex', we expect the quench to lead 
to quadrupolar density oscillations and current reversals, 
similar to what we found for the 'long vortex' in the stripe 
flux superfluid. Inspired by recent work on the superfluid 
Hall effect of atomic bosons [67 , we study the behavior 
of the quadrupole moment Qxy{t) = '^j,xyn{T^t). We 
find that Qxy{t) indeed displays oscillatory sign reversals 
and, as seen in FigjTj^a), the data for various L collapse 
when plotted as Qxy{t)/L^ versus t/L. The t/L scaling 
shows that the oscillations occur due to transport across 
the system length L. A simple scaling argument for an 
edge current induced oscillation shows that Qxy ~ L^, as 
we also see numerically. The numerical observations are 
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thus consistent with the quadrupolar oscihations being 
driven by edge currents. For (j) = — 27r/3, Qxy{t) has the 
opposite sign. Taken together, these observations pro- 
vide strong evidence that the initial state is a nontrivial 
insulator that is incompressible in the bulk and supports 
chiral edge currents. 

In addition to this quadrupolar oscillation arising from 
edge modes, we find that the density in the bulk is also no 
longer constant following the quench. Instead, quench- 
ing Jx leads to x-oriented stripe patterns of the density, 
which originate at the edge and appear to propagate in- 
ward into the bulk. Such a breakdown of the incom- 
pressible quantum Hall state in the bulk can be under- 
stood physically by analogy to the physics of continuum 
Landau levels. A uniform magnetic field in the contin- 
uum can be modelled in the Landau gauge where we set 
Ay = 0^ Ax = By corresponding to a magnetic field Bz. 
This leads to the Hamiltonian 

H^^ = ^/y + ^JP^ + <lByf + Viy) (28) 

where q is the charge, and for simplicity, we assume a 
confining potential V{y) only along the ^/-direction and 
periodic boundary conditions along the x-direction. (In 
cold atom systems, where the gauge field is produced ar- 
tificially, only the product qB is physical and tunable.) In 
the absence of the confining potential, the eigenstates of 
this Hamiltonian take the form ^^^/^(x) = e'^^^^n{y — y^) 
where yk = —k/qB. Here is the nth eigenstate of 
a harmonic oscillator with an energy (n + 1/2)co'l, with 
00 L = qB/ ^myTrix^ and the particle in such a state is 
localized to the vicinity of y = y^. If the confining po- 
tential is varying slowly, with dV/dy <C ool/^b where 
= is the magnetic length, the eigen- 

states remain nearly unaffected by the confining potential 
(up to a small shift of yk) while the energy gets a correc- 
tion V{yk)' Imagine now quenching the dynamics in the 
x-direction by sending the effective rrix oo suddenly. 
This is analogous to sending J^^ ^ on the lattice. The 
Hamiltonian after the quench then takes the simple form 

Ht>o = + V{y), (29) 

which describes a free particle in a potential. Clearly k 
remains a good quantum number. This means that for 
each /c, there are particles initially localized at different 
points yk and described by an initial wavefunction ^n{y— 
yk), which at time t > are free to roll down the valley 
of this potential and delocalize (spread out). At time t = 
0, the particle density is uniform, so there are particles 
localized in the bulk starting at = and going out 
all the way to the two edges at yk ~ ±L/2. For time 
t > 0, this state evolves in time, and the density of the 
resulting state acquires modulations. This leads to stripe 
modulations, with the net density depending on y but not 
on X, with a time dependence governed by a combination 
of the two effects above — particles 'rolling down' the 
potential and spreading of the initial harmonic oscillator 



wavepacket. This picture qualitatively accounts for the 
appearance of stripe modulations of the density in the 
bulk in the quenched state. It also suggests that the 
stripe modulation does not contribute to the quadrupole 
moment dynamics, which is purely an edge current effect, 
consistent with the scaling in our numerical results. 

It is natural to ask how the quench dynamics are mod- 
ified when the boundary conditions are no longer of hard 
wall type, and the atoms are confined to a harmonic trap. 
In this case, we continue to expect bulk stripe-like den- 
sity waves to emerge. However, the density dynamics 
close to the 'edges' that lie parallel to the quenched hop- 
ping direction (i.e., the top and bottom sections of the 
trap) is expected to change. Instead of the density pil- 
ing up at the boundary and reversing, the current will 
continue to fiow up the trap potential before eventually 
reversing. Hence, the frequency of the Qxy oscillation 
will be determined by the trap frequency and the effec- 
tive mass of the fermions, and should match that of the 
bulk stripe density oscillations. Despite this, there will 
remain a definitive signature of edge state currents that 
can be found in the enhanced amplitude of the density 
modulations along the top and bottom edges. (These 
edge states can be thought of as having a finite initial ve- 
locity perpendicular to the quenched hopping direction.) 
This will lead to a shearing of the cloud density and a 
finite Qxy, which is expected to be easily discernible. 

We expect Chern band insulators to exhibit similar 
quadrupolar density oscillations arising from the currents 
at the edge. For quantum spin Hall insulators, with a 
conserved Sz magnetization, we can imagine doing a sim- 
ilar quench experiment and measuring — Q-^y , which 
would exhibit oscillations with a similar scaling. Such 
spin resolved quadrupole measurements rely upon the 
recently demonstrated experimental ability to measure 
spin resolved densities [37] . 



VII. SUMMARY 

Atomic bosons and fermions in the presence of frustra- 
tion or background synthetic gauge fields carry mass cur- 
rents with diverse current patterns or even form gapped 
topological phases with edge currents. We have shown 
that anisotropic quantum quenches can yield a power- 
ful probe of such equilibrium current patterns of atoms 
in an optical lattice by converting them into measurable 
real-space density oscillations. In order to avoid exciting 
particles into the high energy bands of the periodic op- 
tical potential, the quench must be "adiabatic" on time- 
scales comparable to the inverse interband gap, while also 
being "sudden" on time scales governing intraband dy- 
namics. This requirement can be easily fulfilled in ex- 
periments since the tunnel coupling between neighboring 
wells is exponentially suppressed when the lattice depth 
is increased, while the energy separation between bands 
grows with the square-root of the lattice depth Re- 
alizing our proposal for an experimental probe of currents 
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would open up a new avenue to study exotic phases of 
ultracold atomic matter. 



Numerical evaluation of the Gross-Pitaevski 
equation 
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VIII. APPENDIX 



The time-evolution of the initial Hamiltonian's equi- 
librium state after the quench is obtained by numerically 
integrating the respective time-dependent CP equation. 
Specifically, Eqn.[5]is discretized into small time steps Jdt 
that were typically about ~ 10~^. The time evolution of 
the initial state is then determined using a fourth-order 
Runge-Kutta method. In order to ensure convergence, 
this process is repeated many times with increasingly fine 
discretization to confirm that there are negligible differ- 
ences between the solutions. As another check, the to- 
tal energy and particle number are computed to confirm 
that they remain constant throughout the time evolu- 
tion. Eventually, every 1000 steps the density profile and 
other observables are computed using the wave function 
at that instant of time. 



A. Numerical solution of equilibrium GP equation 
for trapped Bose superfluids 



Simplifying the GP equation for quenching of 
the checkerboard flux superfiuid 



For Bose superfiuids in the presence of a trap, the 
ground state is determined numerically using a self- 
consistent minimization of the initial energy functional 
in Eq. |4] recast in terms of the local mean-field density, 

r,r' r 

+C/^nr*:v|/r-|^n2. (30) 

r r 

Here is the condensate wavefunction at lattice site 
r = {x^y) and = |^rP is the average particle den- 
sity at site r. The self-consistent solution of Eqn. 
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IS 



one where the density distribution computed using the 
eigenfunctions of Hmf equals the density distribution 
in the Hamiltonian. In order to find the self-consistent so- 



lution of Eqn. [30] corresponding to the initial pre-quench 
state, we follow a simple iterative procedure. We start 
with a trial density distribution ni^'' corresponding to 
a Thomas-Fermi profile and solve for the single parti- 
cle ground state eigenfunction of Hmf (^r^^ ) • The corre- 
sponding many-body condensate wave function is simply 
given by the normalized single particle solution with a 
multiplicative factor of a/M, explicitly ^r^^ = a/MV^f^^ 
(where M is the number of lattice sites). The condensate 

density is then determined by rir^^ = l^r^-^ p. This is used 
to generate a new trial density distribution via the rela- 
tion ni^"^ = (1 — a)nr"^^ + ani^\ Here, a is strategically 
chosen from (0, 1) to 'throttle' the iterative process in 
order to help maintain convergence and avoid runaway 
solutions. These steps are repeated with the new trial 
density distributions and iterated until the local density 
converges to within 10~^ average variation in the density 
at each site between successive iterations. 



Setting the GP wavefunction to be A{t) + r]rB{t)^ we 
can substitute this into the full time dependent GP equa- 
tion to obtain equations of motion for the complex coeffi- 
cients A{t)andB{t). To obtain the equation for the stag- 
gered density Auab after the quench, it proves simpler 
to define the following variables: 

AuABit) = 2[A''{t)B{t)^A{t)B%t)] (31) 
JC{t) = \A{t)\'-\B{t)\' (32) 
J{t) = i[A%t)B{t)-A{t)B''{t)]. (33) 

Here JCisproportionaltothebondkineticenergyandJ is 
proportional to the bond current. We then obtain the 
equations 

dAuAB 



dt 



dJC 

~dt 
dj^ 
dt 



470/C + 4eo J, (34) 
-ToAnAB - U J Auab, (35) 
-cqAuab + U AuabJC, (36) 



where we have suppressed the time label for clarity. Go- 
ing to second order in time for Auab yields 



= -AX'qAuab 



4U{Jio-JCeo)AnAB. (37) 



d'^AriAB 
dt 

where Aq = 70^ + eo^. To make progress, we resort to 
the following approximation, which is valid at early times 
where we expect well defined oscillations of Auab {t) • We 
replace JC and J7 by their initial values obtained from 
^4(0) and B{0) that correspond to their equilibrium, pre- 
quench, values. Let us call these JCq and jTo- Then we 
find 



d^AriAB 
dt 



4Ag + 4/7(Jo7o - /Coeo) Auab- (38) 
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This yields 



dt 



-Q'^Auab^ where, using the ex- 



plicit values of /Co and j7o, we obtain 
Uno 



|Ao| 



(eoeo +7o7o)) 



(39) 



Since we start at time t = with a uniform superfluid 
having no density modulations, the solution to this takes 
the form An^^(t) = r sin Qt. To find r, we use the initial 
rate of change. 



{dAuAB / dt) 



t=o 



rl] = 4(7o/Co + eoJ7o). 



(40) 



This can be simplified to rQ = 4X((5J/J), where X is the 
magnitude of the initial equilibrium current, which is the 
same on all bonds, and SJ is the amount by which we 
quench the x-bond hopping, leading to the final result 



AriABit) 



(41) 



For a weak quench, SJ <^ J, while a strong quench entails 
setting SJ = — J. The strength of the quench determines 
not only the amplitude of the oscillations, but also their 
frequency We find that this result fits very well the 
early oscillations of Auab obtained by a direct numerical 
solution of the GP equations not only in the continuum 
but also in the central region of the trap for weak as 
well as strong quenches. At later times, the dephasing 
of the oscillations in the trap leads to a decay of the 
Auab arising from spatial variations of Vt via its density 
dependence. 



D. 'Thermal' noise 

By imprinting random phase fiuctuations on the initial 
state, we increase the energy of the system. We expect 
that this noise will lead to an excess energy density that 
will scale as AE ~ ^6>^ax- To check this scaling, we 
computed the ratio of the excess energy for two different 



values of (^^max, choosing S9, 



(1) 



0.5 and Se^^L 



1.0, 



2. If we time- 
evolve this initial state (without making a quench), we 
expect this excess energy to lead to a typical state from 
a 'thermal ensemble' — detailed issues regarding ther- 
malization will be discussed elsewhere. To provide a 
crude estimate of the effective temperature of this 'ther- 
mal state' before the quench, we assume that the domi- 
nant excitations in the system induced by such random 
imprinted phase fiuctuations are the low energy linear 
sound modes. For <C J, the low energy Bogoliubov 
sound mode in the presence of staggered fiux [54 may 
be approximated as hw\^ ~ ck, with the sound speed 
c ~ ^JnU /m* . Here, n is the density, and the inverse 
effective mass is 1/m* = 2Ja^cos(0/4) (where a is the 
lattice constant). Computing the excess energy density 
in the center of the trap, we can estimate the tempera- 
ture of the 'thermal state' as SE = [C(3)/7rc^]T^. For our 



parameters (~ 20 atoms per well, U = 0.2 J and ^ = 7r/2, 
and averaged over states with different initial random- 
ness), we find the temperatures in the two cases to be 
T(i) ^ 3.4J and T^^) ^ 5.2J, which are both signifi- 
cantly smaller than the Berezinskii-Kosterlitz-Thouless 
transition temperature, which can be roughly estimated 
to be Tbkt ^ 7rn/2m* ~ 40 J. This is consistent with 
our assumption that only low energy sound modes are 
excited in the thermal state. 



E. Quench induced density dynamics for the 
staggered fiux state of fermions 

We begin with the staggered fiux Hamiltonian in mo- 
mentum space, 

= ^'(^k<(cos^kr" + sin^kr^)^^, (42) 

k 

where Q = (7r,7r), = (/k^/k+q)' '^^'^ Pauli 
matrices. The prime on the momentum sum implies that 
only momenta in the reduced Brillouin zone are included. 
Here, we have defined l^k = V^k + 7k' cos^k = ^k/^k, 
and sin^k=7k/^k7 with 



^k = —2{Jcoskx 



- Jy COS — COS ky) 



7k 



-2Jy sin — cos ky. 



(43) 
(44) 



This leads to mode energies ±l^k in the initial state. A 
translationally invariant quench of the hopping (say Jx) 
ensures that different momentum pairs (k, k + Q) stay 
decoupled from each other. Nevertheless, this quench 
instantaneously changes ^k ^ s'k, and 7k ^ 7k, so that 
we modify (l^k,6>k) ^ (^k,^k). 

Let us define the initial quasiparticle operators ai^2 
and the final quasiparticle operators Pi^2 via 



and find ^/AE^^y/AE^ = 1.9, which matches closely ^^^^ 

(2) (1) ' 

to the expected value of SOrnL^/SOnii ~ "'^ ^-iw.^ 



sin((9k/2) cos(i9k/2) 
-icos(6>k/2) zsin(l9k/2) 



sin(6'k/2) cos(6'k/2) 
-zcos(6'k/2) zsin(4/2) 






(45) 



. (46) 



Here, the quasiparticle ak,i (<^k,2) of the initial Hamil- 
tonian has energy — l^k (+^k), while the quasiparticle of 
the final Hamiltonian /3k, i (/3k,2) has energy — f^k (+^k)- 
For simplicity, let us assume that we are at a filling of 
less than one fermion per two sites, so that only some of 
the ai quasiparticles are occupied initially, while none of 
the a2 quasiparticle states are occupied (although this is 
easily generalizable to greater fillings). 

We can first transform this into the /3-basis to get the 
dynamics via 



/i[ = sin(^~k/2)e-^^>'*/?^,- 



-cos(^^k/2)e*^''*/?^2 (47) 
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To compute the expectation values, we then need to 
transform back to ai^2 quasiparticles, keeping in mind 
that the ground state at t = has no a2 quasiparticles. 
This means that it suffices to set /3k, i = Q^k,i cos(^k — 

^k)/2 and ^k,2 = c.k,i sin(^k - ^k)/2. 

To compute the density modulation between the two 
sublattices at a subsequent time, AnAsit) = {tia — ^b)^ 
we write 



where M is the number of lattice sites. Using Eq. (47) 
and Eq. (psl), we find 



AriABit) = ^Y^ sin(^k - 4) sin(2(^kt) (50) 



{k}o 



M 



E«/k/k+Q)t + (/,I+Q/k)t) (49) 



where the momentum sum runs only over initially occu- 
pied states in the reduced Brillouin zone. 
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